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NOMENCLATURE 

A = a coefficient used to denote the mapping; function 

2 

as W = - A. Z to represent flow in the corner 
q/ (-1 . \/^. = discharge coefficient 

D = diameter of the circular crest or a parameter 
assumed to describe spillway shape 
g = gravitational acceleration 
H - total head 

h = head over the spillway 
k = boundary roughness scale 
N 3 = number of stream-lines 

n = number of stream-tubes into which the bigger mesh 
near the comer is subdivided 
1 = height of the spillway 
p = pressure at any point 
q = discharge per unit width 
u = velocity component in the x-direction 

V = velocity component in the y-direction 

V = complex velocity 

W = complex potential 

X, y = Cartesian coordinates 
Z - complex variable 

^ = unit weight 

S = angle of inclination of the velocity vector with 
x-direction 



dynajiic viscosity 
density 

surfa.Ge tension 
velocity potential 
strean function 



ABSTRA.CT 


The numerical method of solution of irrotetional flow 
over a spillway, as developed hy Cassidy, v^’as found to he ambigu- 
ous and was not found to converge. 

The method of solution developed hy Cassidy was studied 
in detail, some of the discrepancies were rectified, and a 
workable method of solution has been developed. The method 
developed involves trial and error as well as successive appro- 
ximations. Formulation of algorithms for the correction of 
assumed free surface and assumed value of coefficient of dis- 
charge and refinement of existing method are the main concern 
of this report. 

The coefficient of discharge, free-surface profile, and 
the velocity and the pressure distributions for any given head 
over a given spillway could be obtained by the method that has 
been explained in detail in the thesis. 



CHAPTER 1 


IITRODXJGTIOH 


Spillway performance under varied operating conditions 
is of utmost importance for spillway designers. It is well 
known that when an ogee spillway is operated at heads higher 
than the design head, pressures lower than atmospheric are 
developed on the surface of the spillway and the coefficient 
of discharge increases upto a certain value of the ratio 
of actual head to the design head, heyond which separation 
occurs and the coefficient of discharge decreases rapidly (1). 
This characteristic is taken to advantage by spillway designers. 
Spillway profile is designed for heads less than maximum head 
so that at maximum head a higher coefficient of discharge can 
be obtained. This results in a shorter spillway for a certain 
height of the spillway or a higher spillwa 3 '’ for a certain length 
of the spillway as compared to what would have been obtained by 
taking the maximum head as the design head. The important cri- 
terion of the design, however, is to chock that cavitation does 
not occur due to the development of negative pressure. 

In order to adopt the above mentioned rational design 
procedure, it is essential to investigate the flow- characteris- 
tics for spillways at heads other than design head to determine 
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■the Variation of "the coefficient of discharge and mi n i ngini 
pressure developed on the spillway with the ratio of actual 
head to a reference head (2), 

O The flow over a spillway is of highly contracting nature 
and such flow can be assumed to he irrotational as long as 
there is no separation. Such a flow-problem is well-suited 
to analysis hy the technique of potential theory and the 
results thus obtained would be very near to the actual flow 
over the prototype structure. 

Solution of free-surface was first performed by A,Lauclc 
in 1925 for the case of Sharp Crested weirs of infinite height 
by successive approximation in the con5>lex potential plane (3). 
Using a pair of integral equations derived from the real and 
imaginary parts of Cauchy’s integiral formula, together with 
the integrated form of Bernoulli equation it was possible to 
obtain a new distribution of the free-surface on the basis 
of an assumed distribution. The process was one of successive 
approximation. It was found that the results converged for 
only a particiiLar value of discharge, 

lUV. Southwell solved the problem for a flow-region 
bounded by horizontal floor and free-surface by solving the 
Laplace equation (4). The method required adjustment of the 
assumed free-surface by trial and error and the relaxation 
solution was performed for each trial until the fiee-surface 
satisfied Bemoulli equation. 
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In 1964 Strelkoff published, a paper on solution of 
highly curvilinear ^ro.vity fT-ows (5) but the method suggested ' 
by him is restricted to straight solid boundaries only. He 
solved the flow over sharp crested weir of finite height by 
the use of conformal mapping technique. An algorithm for 
correction of total head vras formuleted by him. 

G-, Watters and R.I. Street solved for the free-surface 
profile over stops and curved hunps (6). They solved the 
integral equation by expressing the complex velocity as in- 
finite series. In neither Strelkoff's method nor that of 
Watter and Street there is any provision of specifying 
independently the shape of the solid boundary. 

Recently a method wa.s developed by J.J. Gajssidy for 
the computation of coefficient of discharge, free-surface 
profile, pressure distribution, etc. for flow, over spillwa.ys 
of finite height (7). A frce-surfaco is assumed initially 
and by an iterative procedure for sa,tisfying the various 
pertinent equations and hounda,ry conditions, the solution is 
obtained. An algorithm v/as obtained to correct the assumed 
value of the coefficient of discharge through the observa- 
tion of behavior of the free-surfa.ee profile from iteration 
to iteration. Ca,ssidy claims this method to he one of 
successive approximations and according to him it is superior 
to other methods involving trial and error procedure. It 
appears that according to Cassidy, irrespective of the 



free-surface assumed initially, the solution converges pro- 
vided the discharge coefficient is correct. It has been 
found from the present study that the eonvergence or diver- 
genee of the solution, obteined by using the method suggested 
by Cassidy, depends on the initially assumed free-surface as 
well as the discharge coefficient. 

The present study was originally intended to apply 
Cassidy’s method for solving the problem of flow over ogee 
spillways with various ratios of design head to the hei^t 
of the spillway. But the available time was spent in supply- 
ing the method only to the case of flow over a spillwo,y with 
circular crest, discovering in this process the drawbacks in 
Cassidy’s method and developing a method of solution. The 
various alternative ways for the computation of the different 
parameters in the procedure suggested by Cassidy have been 
analysed and discussed and a suitable method of solution ha,s 
been developed. 



GHIPTER 2 


MILYSIS OF THE PROBLEM 


2. 1 Physical Analysis : 

If the flow approaching the spillway is ossxuned to he 
suhcritical, i.e. if the Froude Ktunber of the approaching 
flow is less than unity, the control section for the flow 
is at the spillway. Under such conditions the pertinent 
variables that describe the discharge over a two-dimensional 
spillway and the free-surface profile ma,y be grouped as shown 
below: 

q = f^(D, P, h, k, V , p , (f) ...(1) 

^o ^ y » f j » <r) ...(2) 

0,1r,P, V, e = f^^4^5^g^.^(x,y,]>,P,h,k, f,^,0 ( 3 ) 

wherein 0 and f are the two-dimensional velocity potential 
and stream function^ respectively, x and y are the coordi- 
nates of any point, and y^ stand for the x and y coordina;- 
tes of free-surface. V and 9, refer to the magnitude and 
inclination of the velocity vector with reference to x - 
direction, p represents pressure at any point, The para- 
meter D describes the shape of the spillway. P is the height 
of the spillway and h is the head over the spillway crest, 
k is the bounda.ry roughness characteristic, y, P and <S~ 
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are specific weight, density, dynanic viscosity and surface 
tension of the fluid, respectively, g is the gravitational 
acceleration, 

0 and Ip are assumed to be zero at the stagnation point 
B CFig. 1). For the solid boundary -which is also a stream 
line, U is assured to bo zero and for the free~surface Ip is 
equal to the discharge q per unit idLdth. The two equipoten- 
tial lines on the upstream and dovm stream ends of the flow 
region are assumed to form the upstream and do-wnstrean bound- 
aries of flow region under consideration and the stream lines 
are considered to be straight and parallel to the solid bound- 
ary at all points along those equipotential lines. 

X refers to the horizontal distance from the stagnation 
point B and is measured positive do-wnstrean. y denotes the 
vertical distance from the stagnation point B and is measured 
positive upwards. The parameter D is assumed to describe the 
shape of the spillway completely. This type of single-para- 
meter is not obtained for all shapes of spillways. In the 
case of standard spillway this parataotor is the design head 
and in tho case of circular crest spillway this nay be taken 
as the diameter of the circular-crost . In the present analy- 
sis D is assumed to be a length param.otfc.r. 

Using 7T -theorem the folio-wing relationships nay be 
obtained. 




= fg C~J ^ 
^ D D D 


5 , q. . 

CP+h)jc<^7p D) 
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= f f-O ^ JL JL - a -MJ- .. >1 

9 ^ D ’ D * D ^ D ’ /-n ’ 7rT7^rr~-9ir^ ^ 

(P+h)^ (P+Ji)U'7-/f D) 

0 f p V 


■'/^ D- 


?72 ’ Y2g ^ ■’ ^ 

f Z £ k ^ ^ ^ 

^10,11,12,13,14 D’ P’ D» P’ (p+j^)4r 


(P+h) V(r-/f D) 


The discharge and the free-snrface profile are thus 
functions of the relative height of spillway P/D, the relative 
head over the spillway h/P, the relative roughness k/P, a 
Reynolds numher and a Weher number. 


If, on the other hand the Proude number of the approach- 
ing flow is greater than unity and the spillway height is low 
enough to prevent the formation of a surge, the control section 
is upstream from the spillway and the discharge coefficient for 
a given head is automatically determined by the Proude number 
which is an independent variable under this condition. In such 
case, the discharge coefficient is independent of the profile of 
spillway. The analysis reported herein is restricted to that for 
the sub critical approaching flow. 

It has been found that the effect of boundary layer on the 
free-surface profile and the coefficient of discharge can be neglec- 
ted (7). Hence the Reynolds number and the relative roughness may 
be dropped from eqn.(6).In addition; Weber number may be dropped 
if the analysis is — 



restricted to large-scales where the surface tension has 
insignificant role* . H&nce. , 1 equations 4, 5 and 6 are 

simplified to. 


V?ghV2 - ^15 ^ D ’ D ) 

- -F /' _fj2. P h N 
D 16 D » D ’ 1 ^ 

^ p Y 

fzi 3)5/2 ’y^ 1^/2 ’ “5- ’ 




f -S_ _i_ 

17,18,19,20,21 ^ I) » 1 





..(9) 


2.2 Mathematical Analysis 

As the flow is highly contracting and over a spillway 
with a streamlined profile, it maj^ he assumed to he irro- 
tational and so potential theory may he used for the solution 
of the problem {ti) . 

It IS found to he simpler to work in the complex- 
potential plane. When the region hetv/een the solid boundary 
and the free-surface is mapped in the complex-potential plane, 
the physical picture shown in kig. 1 will be transformed to 
a simple rectangular configuration as shown in ?ig. 2. In 
the complex-potential plane the velocity potential 0 is re- 
presented as abscissa and stream function f as the ordinate. 
The complex potential W is represented as 
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¥ = 0 i 


As W is analjrfcic. 


^ _ 
dZ “ 




...( 10 ) 


...( 11 ) 


Z refers to the complex physical plane and is represented 
as 


Z = X + ly . . . ( 12) 

By the definition of 0 and l|l , 



Substituting this in Eqn. 11, 


= 

dZ 



= u - 1 V 


...(15) 


in which u is the horizontal component and w is the vertical 
component of the velocity so that, 

u = V Cos 9 

V = V Sin 6 . . . ( 16) 
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Thus, 


aw 

dZ 


u - 1 V 

V Cos 0 - iV Sin a 


= Y (Gos 9-1 Sin 9) 


= V e 


-10 


...(17) 


Talcing natural logarithm of the q.uantities on both sides 
of Eqn. 17 , 

In (-||) = In V - i0 ...(18) 


Now In (dW/dZ) is anal3rtic, so its real and imaginary 
parts In Y and 9 must satisfy the Cauchy- Riemann equation 
and Laplace equation. 


Hence, hy Laplace equation 

In (V) ^ 8^ In (Y) 
9ljF^ 30^ 

^ = 0 
81|j 2 d02 

and hy Gauchy- Riemann eqn. , 


= 5 (7 6) = 

e In (V) _ a(-e) 

8 ^ 




...(19) 
. ..(20) 


. • * ( 21 ) 
.. .(22) 
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The connection hetvreen complex-potential piano and 
physical plane is formulated hy considering the expressions 
for the stream function and the velocity potential, 


dT^ = 

II '1=^ + ly • 

dy 

...(23) 

d0 = 

H 

dy 

... (24) 


and 



Hence, 

dTjl = - vdx + udy ...(25) 



d0 

= udx -f vdy 

CM 

• 

« 

Solving these 

two 

equations simultaneously. 



dx 

_ ud0 - vdllJ 

2 2 

U 4- T 

...(27) 


dy 

_ vd(2l - ud,$ 

2 2 
u 4- 

... (28) 

Substituting u = 

V Cos 0 and V = V Sin 9, 



dx 

dy 


Oos 9 
V 

Sin e 

— 


d0 - 


Si n 9 
Y 


dljl 


a0 + dijr 


...(29) 


... (30) 



Differentiating the above two eciuat ions with respect 
to 0 and itr a set of four equations is obtained a.s 


dx _ Cos 9 Sin 9 dt 

d^ " V ' ■ Y 

dy _ Sin 9 . Cos 9 dt 

d^ - — ^ 

dx _ Goa 9 Sm 9 

^ “ V If ~ V 

dy _ Sin Q 6.0 , Cos 9 

df - ~T*" If + 


. . ^31a) 


. . . (52a) 


.. .(31b) 
...(32b) 


Expression for distances along the stream line is 


ds . 1 

d^ Y 


. . .(32c) 


finally, for the free-surface the pressure is zero. 
Hence, Bernoulli equation for the free surface is reduced 
to the form given below: 

+ y = H ... (33) 


Thus there are five independent equations out of ten 
equations 19 through 22, 51a, 31b, 32-i, 321, 52c and 33 to 
solve for the five unknowns 9, Y, x, y and H for the known 
values of 0 and f in the 0 - f plane. Only two equations out 
of the four equations 19 though 22 and two equations out of 
the five equations 31a, 31b, 32a, 32b, 32c are independent. 



In order to make the solution as general as possible 
all the ■variables were made dimensionless as described in the 
first part of this cha,pter. How onws-rds all the variables have 
been used to represent their dimensionless form. 

For the numerical solution the flow-region in the com- 
plex-potential plane was divided to form a, grid as shown in 
Fig, 3» Ihe two meshes near the stagnation point were sub- 
divided to form still smaller meshes for the reo.sons mentioned 
leter. 

The following finite difference equation, derived by 
W. G. Biekley (8) was substitiited for equation 20 a.s 

209^ = 4(e^+92+e^+9^) + (e^+egte^+Gg) ...(34) 

In the above equation the subscripts denote the values 
Corresponding to particular node point as shown in Fig. 4. 

Equations 21 and 22 were represented in the integral 
form as, 

A In (V) = A- II i0 ...(35) 

Ain (V) = / ^ ^ dljr ...(36) 

in which A In (V) is the difference in In (V) over the 


interval of integration. 



u 


Eq.uations 31a, 
in integral forms as 

4 X = 

X = 

AJ = 

Aj = 

AS = 


31b, 32a, 32b sDid 32c are represented 


r^2 Sin e 

at 

/z Co^ 

^2 Cos 9 

1 

/s Sin 9 .w 

f -~Y~ d0 

5^1 

X ^ d0 


dijr 


^1 


...(37a) 


...C57b) 


... (38a) 


... (38b) 


...(38c) 


The above equ?,tions ere valid, only a.long Imes of 
const rjit 0 or T|F. 

Mapping function for Singular Point • 

At the stagnation point B a singularity exists. At this 
point the velocity becomes zero and hence the equotions 37 
and 38 for Ax, Ay or As cannot be used at this point. The 
direction of velocity vector is undefined et this point. This 
difficulty was avoided by using the mapping function 

W = ~ A . ...(39) 

near the stagnation point. The coefficient A was evaluated 
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"by equating at point F (Fig. 5) the magnitude of Tolocity 
obtained from equations for corner flov-r and the outside flow, 
Separa.ting the real and imaginary po.rts of equation (39) 

-4" = - + y^ ...(40) 


A 


-2 xy 


. .. (41) 


By solving these tvro equations simultaneously, expressions 
for X and y are obtained as 




Differentiating equation 41 with respect 




Hence expression for 6 is obtained as 
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Using Uqn. 15 ajid differentiating W with respect to Z, 


dT 

= -2AZ = u - iv 
separating the reel part, 
u = - 21x 

suhstitubing the expression for x obtained from Bqn. 4-2 

2A 


■u 


Y = 


Cos 9 Cos 0 


-0 + 


2A 


...(46) 


Solving for A, expression for A is ohtrined as 

2(-04-f 02T'f2y’ 


..(47) 


Y/hile using 3q.n, 34, the sta.gnation point was a. voided 
hy computing the values of 9 for the points adjneont to the 
stagna.tion point by using Eqn. 45. Value of A was not needed 
for this computp.tion. It was computed by Eqn, 47 after the 
velocity for the points ahjacent to the stagnation point is 
known by using equation 35 and 36. To sta 3 rt with the com- 
putation of X and y from the origin, equations 42 and 43 were 
used in place of the general equations 37 and 38, 
















CHilTSR 3 


NUMERiaiL SOLUTION 
Nunorical Method; 

The naln am of the problaa is to solve the equations 
along with pertinent boundary conditions developec* in 
chapter 2 for obtaining tho flow net and the coefficient 
of discharge for irrotational flow over a spillway. 

Various alternative procedures to satisfy the 
Various equations have been tried. The one found to 
be successful has been described in this Chapter. The 
behaviour of the solutions observed in the rest of the 
procedures studied has been discussed in chapter 4. 

ils stated earlier all the ccnputations have been 
made in the Cvmplox-petential plane. For the given 
spillway shape and head cvt,r the spill^'/ay a rough free- 
surfacG was obtained by drawing a flow net for assuned 
free-surfacc profile and - 
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correcting th.e free-surface on the hasis of the flow-net drawn. 
This was done to promote convergence. 

On the basis of this rough flow-net obtained, values of 
0 were assumed for all the nodal points on the boundaries in 
the 0-ljr grid as shown in figure 3. By using equations derived 
in previous chapter it was then possible to compute all the other 
variables 0, V, x, y etc. for all the nodes in the congplex-poten- 
tial grid and thus it was possible to obtain a new distribution 
of 9 for the boundary. This procedure was repeated for a few 
iterations (generally two or three) to observe the behavior of 
the free-surface. By observing the behavior of the free-surface 
obtained for a few different values of C^, covering almost the 
whole range it was possible to make necessary corrections in 
the initial assumption of distribution of 0 of free surface, if 
the same was in large error, later, when the assumed free-surface 
was reasonably corrected, the procedure was repeated for different 
coefficients of discharge and the particular value of discharge 
coefficient for which the change of x and y coordinates of the 
free surface after a certain number of iterations is less than 
1% of the values before tbe iteration was considered to be the 
right value. 

The free surface obtained after each iteration diverged 
upward or downward depending upon the initial assumption of 
free surface as well as the value of coefficient of discharge. 
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It Converged to the right free-surface for a particular value of 
that too when the initial assumption was not very 

wrong. 


The algorithms obtained for correcting the value and 
the initially assumed free-surface are discussed later in this 
chapter. The procedure adopted for correcting the free-surface 
profile on the basis of the assumed or old distribution of 9 
listed below: 

1. The value of discharge q was coii 5 )uted for the assumed 

value of G^, using equation, 


q 



^d 


h 


3/2 


wherein all the variables should be in thoir dimensionless form. 


2. The flow- region was divided into suitable number of stream 
tubes depending upon the required accuracy and available computer 
time. The difference A m the values of l|l for two adjacent 
stream lines is given as 

Af = q/(Ng - 1) 

where N is the number of total stream-lines chosen for computa- 
s 

tion. 

for using equation 34-, it was required to assume square 
meshes in the 0-l]r grid and hence the difference A0 between 
two adoacent equipotential lines was taken as equal to (fig. 3). 
The number of equipotential lines was chosen on the basis of the 
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roTigli flow-net drawn and by assuming the positions of the up- 
stream and downstream boundaries enough to assume the eq.ui- 
potential lines to be straight lines on the tz/o boundaries. 

3. On the basis of the rough flov'~net dravm, distribution 
of 0 was assumed for the nodal points on the boundaries, Ibr 
the upstream and downstream boundaries tliL^se were assumed to be 
constant and eq.ual to zero and slope of the spillway at the 
downstream boundary respectively, 

4. 9- field for the interior points in the complez potential 
grid was computed by successively applying equation 34 to each 
anterior point hy the relaxation procedure, until the v^^lues 
obtained after two consecutive iterations differed within the 
permissible limit. 

5. Velocity distribution along the free-surfoce v/as computed 
such that it satisfied both equations 33 and 38h simultaneously. 
The computation was started from the upstream end where the 
velocity distribution was assumed to he uniform, for the succes- 
siTG points, firstly a rough value of y was assumed and by using 
eqn. 33, V w^^s computed. This was substituted in equation 38h 
to g('t a mortr- correct value of y. Iterations for y and V were 
done till the two successive values were -"ithin specified limit. 
This was repeated for all successive points upto the domistream 


end. 
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6, The velocity distribution was known at this stage for the 
free-surface as computed in step 5 and along the upstream boun- 
dary by assuming uniform flow-distribution along that boundary. 
Hence, velocity at any interior point and along the solid boun- 
dary was computed by integrating along both the equipotential- 
line and stream-line separately by using equations 35 aud 36 
and then taking the mean of the two values of In Y, 

7. X and y coordinates for all the nodes were computed by 
Using equations 37a, 37b, 38a, 38b and 38c. Computation was 
sta.rted from the stagnation point, the origin. To start v/ith, 
equations 42 and 43 were used to avoid the singularity at the 
origin. 

Using equation 38c the values of s(!^, ijJ) were computed 
for each point along the spillvray profile. It was assumed to 
be zero at the origin and v/as taken positive downstream. x(0,l|J) 
and y(0, l{f) were computed along the spillway profile by using 
their functional relationship with s(0, I^) for the given profile 
of the spillway. With the help of eqn. 37b, x(0, 1 ^) was com- 
puted along the bed upto the upstream boundary. y(0,tf'> along 
the upstream boundary was known by assuming uiiform flow. 

x(0, ijl) and y(0, Ijl) for the other interior points and the 
free surface was computed by integrating along the streamlines 
starting from the upstream boundary and also along the equi- 
potential lines starting from the solid boundary and then taking 
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■fclie arithme'fcic mean of flae "two values "thus obfamed. 

8. The "behavior of the free- surface thus obtained -was observed 

and another value of was assumed accordingly. At this stage 
it was also decided whether a correction in the assumed 0- values 
for the free-surface was required or not. If it was found to be 
essential, the correction in the G-distribution for the free- 
surface was done suitabl;^ and the procedure was cycled back to 
step 5. The algorithm for obtaining the rough free-surface to 
use as an initial assumption for correcting the value of as 
well as obtaining a more refined free-surface by the process of 
successive approximation is described below. 

If the negative values of 9 assumed for the free-surface 
for any portion were too high, the free-surface in that region 
obtained after each iteration progressively diverged downward, 
i«e, » 1^) free surface decreased after each itera- 

tion for the whole range of values of 0^ (Pig, 6). On the other 
hand if it was too low, the free-surface obtained after each 
iteration diverged upward, i.e. y^C^j ”^1^) for the free-surface 
increased progressively after each iteration for the whole range 
of values of C^, On this basis it was possible to correct the 
assumed B-distribution by making suitable adnustments. 

After the approximately correct free-surface was obtained 
by the method mentioned above, the divergence in different direc- 
tions, upwards or downwards was observed for the free-surface for 
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different values of G^. If the value of was too high, the 
free- surface diverged upwards and on the other hand if it was 
too low, the free-surface diverged downvrards, as shovm m Pig. 7. 
Only for a particular value of the coefficient of discharge the 
free-surface converged to the correct one provided the rough 
corrections in the O-distrihirt ion are made initially to remove 
the divergence due to wrong initial assumption. 

All the computations were done on the IBM 7044 computer. 

The large storage required for the In (Y) and 0 arrays zeliminates 
the use of small computers. 

For all integrations, Simpson's rule and half- Simp son 's 
rule (Bef. 9) were used. 

Solution has Been obtained for a circular crested spillway 
for head equal to 2.5 D, wherein B is the diameter of the circular 
jrest. The free-surfacc and flownet obtained are shovm in Pig. 8 
nod the pressure distribution obtained along the spillway surface 
.s shown in Pig. 9. Twenty one - ' streamlines and ninety 

equipot ent lal lines v^ere chosen for computation. 
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DISCUSSION OF RESULTS CONCLUSIONS 
'I* 1 Discu ssicn i 

Cassidy has shewn that the divergence cf the free- 
surfece depends only -upon the value cf cc.efficicnt cf discharg- 
But the present study has shewn that the divergence cf the frei 
surface, as reported by Cassidy, depends both upon the initial 
assumed free surface as well as the value cf C^. The divergeni 
was found tr be dependent cn the value cf alone only for a 
few initially assumed free-surfaccs and such a free surface 
cculd be obtained after a few trials. Tlius the nethed has net 
been found tc bo a simple nethed rf successive approximations 
but involving trial and error procedure also . 

In the procedure suggested by Cassidy, the velocity dis- 
tribution along the free-surf\cc is obtained directly frern the 
values cf y obtained at th. end of the previous iteratim by 
using the Bernoulli equation. So, if a free-surface higher tha 
the correct frec-surfaco is assuriOd, a low velocity distributir 
throughout the field is obtained for a particular value of co- 
efficient cf discharge. This results in higher va.lues C'f y for 
the free surface if the integration f - r y is done almg the 
cquipotential line starting fr a the solid boundary. Thus if 



the initially assumed free-surface is higher than the correct 
free-surface , it becomes still higher and similarly if the 
initially assumed free-surface is lower than the correct free- 
surface it becomes still lower after each iteration, even if 
the correct value is assumed for the coefficient of discharge. 

On the other hand, if the integration for y is started from 
the upstream boundary, the behavior of free surface can never 
be in agreement with what Cassidy has shown. In this case, no 
algorithm for correcting either the initially assumed free-sur- 
face or the value of is obtained. 

The behavior of the free-surface profile in the region 
upstream from the spillway obtained after each iteration using 
this method becomes erratic particularly due to the fact that 
even slightest error in the initially assumed ordinates of the 
free-surface results in error in the velocity distribution along 
the free-surface upstream from the spillway which is computed 
directly using Bernoulli equation. ^Vhen y 13 computed for the 
free-surface by integrating along the eqnipotential line com- 
mencing from the solid boundary the error in y is magnified in 
the region upstream from the spillway as the velocity head is 
very low in that region as compared to the head. 

Bor the solid boundary, in Cassidy's method, x(0, f) is 
computed hy using Eqn. 37 and y and 9 are then computed by 
using the equation for the profile of the spillway. This has 
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the great draw back that the mmber of eq_uipotential lines join- 
ing the vertical face of spillway remains same after each itera- 
tion as assumed previously because always remains zero for 
vertical face and is not corrected after any iteration if the 
initial assumption is wrong. 

Cassidy assumes a mapping function for the corner flov? 
independent of the exterior flow, which does not seem reasonable. 
But in the present case the mapping function for the corner flow 
has been obtained from the condition that along a reasonably 
chosen boundary enclosing the corner, the mapping function and 
tho exterior flow eq.uations should give the same results. The 
accuracy of such an assumption may be verified by bringing the 
boundary closer to the corner and observing that there is no 
change in the results. 

In the method developed and presented in the previous 
chapter, all these drav/backs have been removed. The behavior 
of free-surface profile remains more stable and the behavior 
IS not erratic so that the algorithm for the correction of the 
value of may be used easily. 

4, 2 Conclusions : 

In the present study, a method has been developed for 
tho numerical solution of irrotational fLov' over spillways of 
finite height. A full study of the method suggested for the 
same by Cassidy has been conducted and the drawbacks found in 
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that method have been rectified. The method of solution has 
been found to he one involving trial and error. 

The chief advantage of the method developed herein lies 
in the fact that there is no erratic behavior of free-curface 
as seen in the method suggested by Cassidy, if there is any 
error in the initially assumed free-surface. The solutions 
obtained after each iteration remain more stable. 

The greatest disadvantage of the method is the large 
computer time required for the relaxation of 9 by using equation 
34. In the example presented herein, 26 minutes of operational 
time was tahen on IBM 704-4 computer for 3 iterations. The time 
fcnlcen depends very much on the number of stream-lines chosen. 

The following alternate methods of approach are suggested 
for future investigations: 

1. By using equation 19 instead of equations 21 and 22 for 
the computation of velocity and using the Gauchy-Rc-imcnn equa- 
tion for 9 instead of the Laplace equation the method may be re- 
fined and the solution may become more stable (10). 

2. Integration for x and y may be carried along the stream 
linos starting from the downstream boundary also. This may 
lead to some other divergence-pattern. 
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3 The effect of the location of upstream and downstream 
boundary on the results^ particularly on the minimum pressure 
on the surface of the spillway should be investigated. 
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LISTING OF COMIUTim TROGRilWE 



DESCRIITION OF V/JlimES 


A- 2 


CD 

H 

HD 

I 


ITEST 


J 


MAXIT 


NF 

FS 

NFEMID 

KFC 


= coefficient of discharge 
= head ever the spillway crest 
= paraneter D as described in cha.pter 2 
= Subscript used tr denote equipctential lines 
Cl = 1 refers to the upstrean boundary and 
I = NF refers to the downstrean boundary) 

= an index 

(data f-r ITEST is taioen as positive vAien 
integration for y is carried only along the 
equipotential line and ITEST is taken as 
negative when integration is carried through 
the tw: different paths as described in 
chapter 3*) 

= subscript used to denote strean lines. 

(J = 1 refers tc the free surface and J = NS 
refers to the solid boundary) 

= naxinua nunber of iterations for the re- 
laxation of © 

= total nunber cf equipctential lines. 

= nunber cf strean lines 
= value -^f I for = 0 line 

= nunber cf equipotential lines corresponding 
tc the finer grid near the corner 
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NSC 

NFCMID 


S 


THI TA 

V 
W 
X 

Y 


= nuntier cf streaii lines ccrre spending to 
the finer grid near the corner 
= Value cf I fer 0 = 0 line ccrrespondlng 
tc the finer grid near the corner 
= dinensicnless pressure 
= Distance alo ng the surface cf spillway 
startin^. fr-n tne origin 
= angle rf inclinatic'n of velocity v^ 3 ctcr 
with x-directi'.>n 
= nagnitude rf velocity v ctrr 
= height of the spillway 
= X. Co -ordinate 
= y CO — ~ rdin at e 



LISTING OF P'?OCRANME 


$16 JOB 

SIBFTC MAIN 

dimension TH I TA( 140,21 ) ,STHITA( 100,21) , VL06{ 100 » 21 ) ,V( 140,21 ) ,X 
li^O,21) ,Y( 140,21) ,DuMX (140,2) ,DuMY( 14^,2) ,TH I T As { 140 , 2 ) ,YS( 140 ) ,5 
140) ,THITAC(21,11 ) ,VC(l 1 ) ,VLOGC(ll) ,XC < 1 1 ) , yC { 1 1 ) , STHI TC ( 21 , 1 1 > 
dimension p{140) 

read l,NF,NS,NFlMID,MAXlT,f IT,ITEST 

pEAD i,nsc,mfc,mfcmid 

COMVOM THITA,X,Y,NS,NF 
COM’‘’ON THITAM 
COMMON NFIMID 
COMMON HD 
CALL FLUN(50n) 

CALL FL0V(5O0) 

MCHPNT=1 
READ 2,W,H,HD 

C w=W/HD,H=H/HD,HT=HT/HD 

HEAD 3,CHECKT»CHECKX,CHECi<Y,C:D,CDI »CHECkv 
3 FOpNAK 8F10.B ) 

THITAM=~ATAM(1.)^^4./3. 

PI=ATAN ( 1. )*4. 

c thitam=thita at D/s end 

WS1=NS-1 

mF1=MF-1 

fimid=nfimid 

FNS=NS 

FNSl =NSl 

FNF=NF 

FNSC=NSC 

FWFC=NFC 

rC.MID=NFCMiO 

FNSC1=MSC-1 

MFC1=NFC-1 

NSC1=NSC-1 

C MAPPING FUNCTION W=A^fZ**2 

DO 111 = 1, NF 
DUMX ( I ,1)=0. 

DUMX ( 1 »2)=0, 
nuMYd ,l)=o. 

DUMY { I ,2)=o, 

DO 11J=1,NS 
STHITA( I , J)=0. 

THITA( I ,j)=0. 

11 CONTINUE 

00 413 1=1, NFC 



n rs 


DO 413 J=1,MSC 
THITAC( I»J)=n. 

^13 STHITCI I,J)=o» 

^'^LL IMPUT ( Y,THITA,£,mF,MI‘'I1 ,n'P 1 ,H »t«',THITAM ,NFIMID,NS ) 

BOUND (X,Y,THITA , S ,NF I M I D , MF , NS , THI TAm ) 

PRINT 7 ♦ (X ( I ,NS ) » Y( I »NS) , 1=1 ,NF) 

K = 0 

DO 12 J=1 , MS, MSI 
K=K+1 

DO 12 1=1, NF 

12 THITASI I ,K) = THITA ( I ,J) 

DO 13 1=1, NF 

13 YS(I)=Y(I,1) 

320 CONTINUE 

IFIRST=1 
K = 0 

DO 14 J=1,MS,NS1 
K = K+1 

DO 14 1=1, NF 

14 THITA( I ,J)=THITAS( I ,K) 

DO 151=1, NF 

15 Y(I,1)=YS(I) 

3no CONTINUE 

Q = (2,/3, )*CD«-{H^f-"fl.5) 

HT=H+W+-( 0/ (H+W ) ) 

C 0 IS IN DIMFNSTONLFSS FORM 

DP1=Q/{FNS-1. ) 

DP2=DP1*2, 

DN1=-DP1 
[)^yl2 = -DP2 

dpic=ddi/fmsci 

DP2C=DP2/FNSC1 
DM1C=DM1/FMSC1 
nM2C=DM2/FMSCl 
PRINT?, DPI ,DM2,DP1C,DM2C 
IP1=NFIMID+1 
IMl=NFIMID-l 

RELAXATION METHOD FOR THITAS 
THITA FOR CORNER 
THITAC(NFCMID,NSC-1)=PI/4. 

SINE=-1. 

ICM=NFCMID-1 
ICP=NFCMID+1 
DO 54 I=TCM,ICP,2 
FI=SINE*DP1C 



A*€ 

THITA(I'FIMID-1,NS) = 0o 
THITA(mFIMID+1»WS)=PI/2« 

281 CONTINUE 

DO 330 ITM=1,NIT 
DO 70 IT=1,MAXIT 

THITA{MFlMlD-l,NSl) = (TriITA(NFlMlD,NSl )+THITA( NFIivilD-2»NSl )+TriIT 

IF IM ID-1 ,NS)+THITAINFIMID-1,NS-2 ) )/4. 
THITAC11,11=THITA(NFIMID-1»NS1 ) 

THlTA(NFlMID+lsNSl) = (THITA(NFlMlD+2,NSl)+ThlTA(NFI.MlD+l,NS-2)+Tf- 

lAlNFlMlDsNSl )+THITA(NFIMID+1 ,NS) ) /A. 

THITAC(NFC,1 )=tHITA (NFIMID+1 ,NS1) 

THITA(nFImID,NS1 )={THITA{NF ImID+1 »NS1 >+THITA(NFIMID,NS-2)+THITA 

1IMID-1»NS1 )-^THITAC(NrCMlD,2)-i^FNSCl)/ (3. + FNSC1 ) 

THITAC(NFC.mID,1)=THITA(NFIf,ID,NS1} 

THITAC(1,NSC}=THITA(NFIMID*1,NG) 

THITAC{NFC,NSC)=THITA(NFIMID+1 ,hS) 

DO 51 I=2»NrCMID 

THITAC( I , 1)=THITAC( 1-1,1 ]+(THlTAC(NFCMlD»li-THITAC( 1,1 ) >/FnSCi 

51 THITAC( I ,MSC)=THITAC(1,NSC) 

ii=nfcmid+i 

DO 52 1=11, NFC ^ 

THITACc I , l)=TriITAC(I-l,l) + ( THI TAl ( NFC , 1 ) -THI T AC ( NFCM I D , H )/Fh SCI 

52 THITAC( I ,NSC}=THITAC(NFC,N5C) 

DO 53 J=2,NSC ^ ^ 

THITAC(1,J)=THITAC(1,>»1) + < THITAC{1,NSC}-THITAC( 1,1 ) )/rMSCl 

53 THITAC (NFC,J)=THITAC(NFC,J-1 ) + (THITAC(NFC,NSC)-THITAC( NFC,1 ) )/r\ 
11 


62 

63 

64 
61 


60 

50 


54 


DO 50 J=2,MS1 

DO 60 1=2, NFl 

IF(J-MS1)61,62,61 

IF( I-ir'll)63»6C,63 

IF(I-NFIMID)64,60,64 

IF(I-IP1)61,50,61 

CONTINUE 

SUM1 = THITA(IsJ-1)+THITA(I,J+-1)+THITa{I“*1»J)+THITA(I-*-1,J) 

SUM2=THITA( I+1,J+1)+tHITAI I+ l, J-1 )+ThITA( I-1,J+1)+THITA( 

THITAt I ,J ) = (4o^SUMl+SU,T2)/r0v 

CONTINUE 

CONTINUE 

SI=DP1C 

EX=*FI + ( F ) **0«5 


DYDX=SI /EX 

TillTACd ,NSC-1)=ATAN(DyDX) 


SINE=-SINE 

CONTINUE 


) 



462 

463 

464 
461 


i^OO 

70 

110 


D0450 J = 2,/\jSCl 
DC ^60 I=2,NFC1 
IF( J-USC1)461,462,461 
IF( I-iCM) 463 ,460,463 
IF( I“MFCMID)464, 460,464 
IF(I^ICP)461, 460.461 
CONTl.NiUE 

SUM2 = THITAC( ^ » J+I ^+THItAC( I-l ,j;+tHITAC( I + i ,j) 

CJM2 ThITAC(I+i,j+l)+XHlTAC(I+i,j*i)+THITAC(I-l,J+l)+THlTAC(I-l 

Join Nli E ’ ^ • 

CONTINUE 
DO 80 1=2, NFI 
DO 8o J = 2,,NS1 

^ I , J) > /THITA( I . J) 

tRRT=ABS( DIFFT ) 

IFJEPRT-CHECKT)80,90.90 

continue 

DO 480 I=2,NFC1 
DO 480 J=2,NSC1 

^ >J)-THITAC( I ,J) J/THITAC( I ,J) 

ERRT=A8S( DIFFT) 

I F ( ERRT-CHECKT ) 48 0,490 ,490 

CONTINUE 

GO TO 110 

CONTINUE 

CONTINUE 

DO loo 1=1, NF 

DO loo J=1,NS 

oTHlTA( I , J)=tHITA( I ,J) 

DO 500 1 = 1, NFC 

DO 500J=1,NSC 

STHITC( I , J)=THITAC( I,J) 

CONTINUE 

CONTINUE 

COMPUTATION OF LOG OF V 

'^^LL FrISURI Y,THItA.V,Ht.Dpi,CHECKV,K!F ) 

J = 1 


DO 120 1=1, NF 
D2o=HT«Y ( I , J ) 
IF(D201290,111,111 
111 V ( I , J) =D2C^^^0»5 

VLOGd , J) =AL0G{V C I , J) ) 
]2o CONTINUE 
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UVEL=Q/ (H+W) 

1=1 

DO 121 J = l,NiS 
V( I ,J)=UVEL 

121 VLOG( I ,J) =ALOG(V( I , J) ) 

call VElFI { THITA,VLOG,,'|F,nS,NFIMTD»mF1 ,NS1 ,DP1 ,DP2 »DM1 9DM2»PI ) 
VLOGC( 1 ) =VL0G(NFIMID,N51 ) 

I=NFCMID 
J = 2 

D1={THITAC(I+1,j-1)»THITAC( I-l^J-l) j/D^zC 
D3=(THITAC(I + 1,j+1)-.THITAC( I-l.J+l) )/DP2C 
D2=(THIT-aC(I+] » J )-THITAC( I-*1 » J ) ) /DP2C 
DVLnGC= ( 5 o^D1+8o*D2-D3 ) C/12o 

VLOGC ( j ) =VL0GC ( J“1 ) +DVLOGC 
DO 510 J=3,NSC 

D1 = «THITAC( I+l , J-2)-#THITAC( I-l.J-2) )/0^2C 
D2= (THITAC{ I+l , J-1 }-.THITAC( I»1 .J-1 ) ) /DP2C 
D3=(THITAC( I+1,J)«»THITAC(I-1»J)) /DP2C 
DVLOGC= (Dl+4o'>f-D2 + D3 )*DMlC/3® 

510 VlOGC( J)=VLOGC( J-2)+DVL0GC 
DO 1A5 1 = 1, 

DO 1A5 J=1,MS 
145 V( I, J)=EXP(VLOG( I ,J) ) 

DO 545 1=1, NSC 
545 V'" ( I ) = EXP (VL0GC{ I ) ) 

A=(VC(NSCl)r-*2)/(4v,*DPlC) 

PRINT7,A 

C COMPUTATION OF X AND Y NEA'=' CORNER 

XC(NSCl)=-{DPlC/( 2 o*A) i^*0.5 
YC(NSC1)=»XC(NSC1 ) 

I=NFCMID 

J=NSC-2 

D1=-SIN (THITAC ( I , J+] ) >/VC( J+1) 

D3=-SIN (THITAC ( I , J-1 ) ) /VC( J-1 ) 

D2=-5IN(THITAC { I , J) )/VC ( J5 
DX=( 5o*D1+8u^D2-D3)*DP1C/12o 
XC( J )=XC( J+1)+DX 
D1=CoS( THITAC ( I , J+ 1 ) ) /VC ( J+l ) 

D3=C0S(THITAC( I»J-1) )/VC{J-1) 

D2=C0S( THITAC ( I,J))/VC(J) 

DV=( 5.*D1 + 8<,*D2-D3)-^DPj C/l?o 
YC{J)=YC( J+1)+DY 
DO 550 K=3,NSC1 
J=NSC-K 

Dl=-SIN { THITAC (NFCM ID, J ) ) /VC( J) 
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D2=-SIN ( THITAC ( NFCMID* J+1 ) ) /VC ( J+1 ) 

D3=~SIN(THITAC( I » J+2) )/VC( J+2) 

DX= ( 01+A, *02+03 ) ->f DPI C/3 . 

XC{ J)=XC( J+2 )+DX 
03=C0S{TI!ITAC(KFCMI0,J) )/VC( J) 

OA=COS( THITAC (NFCMID* J+1 ) ) /VC< J+1) 

05=C0S(THITAC( I»J+2) ) /VC (J+2) 

DY=( D3+A , ^DA+05 ) +DP1C/3 , 

YC{J)=YC(J+2)+DY 
550 continue 

X(NF1mid,NS1)=XC( 1) 

Y(NFIMIO.NSl)=YC<i) 

X {NFIMID+i,NS)=0. 

Y (NFlMTD-l,NS)=o. 

PRINT 7,X(MFr'ID,NSl) ,Y(NFIHID,NS1) 

PRINT 7, (XC( J) ,J = 1,MSC) 

PRINT 7, (YC( J) ,J = 1,NSC) 

PRINT 7, (VC( J) ,J = 1,MSC) 

call XYC0MP(X»Y»THITA*V,S*NF»NS»NF1»NS1 tNF IMID,DP1 ,DMl,Df 

lDM2,iriRST,ITEFT,H,W) 
print 7,(X(I,NS) »Y(I»N<;),I = 1,NF) 
call B0UND(X ,Y»THITA,S»mFIMI0,mF,NS»TH1TAM) 

K = 0 

DO 2oo J=l,NS*NSl 
K=K+l 

DO 2rio 1 = 1, NF 
IST1=I 

OIFFX=(DUMX{ I »K)-X( I»J) )/HT 
DIFFY=(DUMY{ I »K)--Y(I*J) )/HT 
ERRX=AB5(DIFFX) 

ERRY=ABS(DIFFY) 

IF(FRPX-CHFCKX) 205,210,210 
?05 IFCEPPY-CHECXY )200»210,210 
20(3 CONTINUE 
GO TO 310 
210 CONTINUE 

DO 230 I=1*NF 
DUMX{I,1)=X{I»1) 

DUMYd ,1)=Y( I»l) 

DUMX(I»2)=X(I»NS) 

230 DUMY(I»2)=Y( I»NS) 

C COMPUTATION OF NEW THITAS 

J=1 

DO 241 1=2, NFl 

240 THITA( I , J )=ATAN( (Y(H-1»J)-Y(I-1,J))/(X(I + 1.J)-X(I-1»J))) 
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DO 256 J=1,NS1 
25n THITA{ I ,J)=0. 

256 CONTINUE 
I=NF 

DO 266 J=l,NSl 
26o THITAl I , J)=THITAM 
266 CONTINUE 

PPI''>T 6,IT,IST1*CD 

PRINT 7,(X(I,1),Y(I,1)»I=1»NF5 

PRINT 7 , (X( I ,NS) ,Y{ I ,NS) ,1 = 1 »^'F ) 

PRINT 7 , (THITA ( I ,1) ,1 = 1 sNF) 

PRINT 7, (THITAd 5NS)»I = 1*NF) 

330 CONTINUE 

PRINT 7,( (X(I ,J) 5Y(I»JI»I = l9r'F) ,J = 1,NS) 


J=NS 


PO 340 I=IP1,NF 

p( I )=HT-V( I , J)-»«-2-Y{ I ,J) 

PRINT 6,1 .J,X( I » J) »Y(I »V( I »J) ,P( n 
34n CONTINUE 

DO 350 I=IP1,NF 
VD=DP2/(S{ I+l )-S{I-l) ) 

ER= (V( I , J)~VD ) /V ( I ,J) 

PRINT 6 ,I ,J,V{ I , J) ,VD,ER 


35o CONTINUE 
J = 1 

DO 360 1=1, NF 

P( I )=HT-V( I,J)**2-Y(I »J) 


PRINT 7,P(I) 
360 CONTINUE 

call 


XYCOMP(X,Y,THItA,V,S,NF ,ns 


,NF1 ,NSl »NFIMID,D pi ,DM1 ,DP 


inM2,lFIRST,ITEST,H,W) 
GO TO 281 
280 CONTINUE 
290 CD=CD+CDI 
GO TO 320 
1 F0R,'1AT(2nl3) 

2 FCR!4AT(8F10.5> 

6 FORi'^AT ( 2 I 5,4E15.7 ) 

7 F ormAT ( lOE] 3 ,5 ) 

8 coRriAT(8Fl5.6) 

9 format { 6F13 * 5 ) 

6on FORMAT! 10F8. 5) 

31n CONTINUE 
STOP 


END 
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SIBFTC VELFI 

^subroutine VElFI ( ThiT a, VLOG, N F, N S,NFIMI D, NFI ,NS 1 ,DP1 ,DP2,DM1 , 

iPlsNFlMjo+i 
DO Ao'l 1=2, MFl 
DO 4n0 J=2,MS 
IF( J-2 )20,10,20 

D3-{THITA( I + Uj+l )-.THITA(I^1 ,j + i ) )/Dp2 

D2={THITA(I+i,j)_thitA( I- l,j) )/5 p2 

DVLOG- ( 5.*D1+8.*D2-D3 )*DM1/ 12, 

VLOG { I j J ) =VLOG ( I » J-1 ) +DVLOC 
13 IF(I^2)11,12,11 

12 D1=0. 
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D2=-(THITA(I »J+1 )-THITA( )/DM2 
D3=-(THITA{ I+1,j+i)«THITA( I+ l, j-i) ) /DM2 

31 n=(5.*Dl+8,*D2-D3)*DPl/12. 

VLOG( I »J)={VLOG( I ,J)+VLOG( I-l» J)+D) /2. 

GO TO 4oo 

11 D1=-(THITA( I-2,J+1)-THITA( 1-2 ♦J-D )/DM2 
D2=-(THITA( I-1,j+i)-.tHITA( I-l , J-l ) ) /Dl'i2 
D3=-(THrTA( I ,J+1)-THITA(I »J- 1 ) )/Dm2 

32 D= (D1+4.*D2+D3 )*DPl/3, 

VL0G{ I ,J)=(VLOG( I ,J)+VL0G( I-2»J)+D) /2. 
GO TO 4no 

2n IF ( J-MS)60,7r,6n 

6 0 D1=(THITA{I + i,j-2)-THITA(I,-1,J-2) )/DP2 
D2=(THITA{I+1,j-i )-THITA( l-l, J-1 ) ) /DP2 
D3=(THITA( I+1 ,j)-THITA( I- l ,J) )/Dp2 
GO TO 300 

7o IF( I-NFIMID)80»90,8n 

9o GO TO 4on 

80 IF ( I-IMl ) 100,110. 100 

110 D1=(THITA( I+l , j-2 )-THITA( I-l ,J-2 ) ) /DP2 
D2=(THITA( I+l, J-l )-THITA( l-l,j-l ) ) /DP2 
D3 = 0, 


GO TO 3n0 

loo IF( I-IPl )6n,110.6n 
300 DVL0G=(Dl+4.^tD2+D3)-<fDMl/?. 
VLOG { I , J ) =VL0G ( I » j- 2 ) +DVLOG 
IF( J-NSl ) 14,15,14 
If IF(I-NFIMID)17,16,17 
16 THITA(I,J+l)=PI/4. 
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GO TO 11 

17 IF( I-IPl) 18,19,18 

19 THITA( I-l ,j+l)=pI/4. 

GO TO 11 

18 IF( 1-101-1)13,33,13 

33 THITA(I-2,j+1)=pi/4. 

GO TO 11 

14 IF(J-NS) 13,21,13 

21 IF(I-NFIMID)22,4no,22 

22 IF( I-!P1)23,400,23 

23 IF( I-IP1-1)25,26,25 

25 IF{I-2)27,26,27 

26 D1=-{THITA( I-.1,J)-THITA(I-1,J-1) )/D'U 
D3=-(THITA( I + 1,J)-.THITA(I + 1,J-1) )/DH1 
D2=-(THITA(I ,j)-THITA(I ,J-1) )/Dm1 

GO TO 31 

27 D1=-(THITA(I-2,J)-THITA(I-2,J-1) )/D1'11 
D2=-(THITA(I-1,J)-THITA(I-1,J-1) )/DM1 
D3=-{THITA( I,J)-THITA( I ,J-l) )/D m1 

GO TO 32 
4on CONTINUE 
I =NF 

DO 5on J=2,NS 
5no VLOG( I , J)=VLOG{ I ,1) 
return 
PND 





SIBFTC BOUND 

subroutine BOUND(X,Y,THITA,S»NFImID,NF,NS,THITAM) 
dimension X( 140,21) 5Y(140»21) »THITA{ 140, 21 ) ,S(140) 
J=NS 

PI =ATAM ( 1 * )*4, 

ST1=2.5 

ST2=ST1+PI*5./12. 

XT2 = 0.5*( l«-COS{PI*5./6. ) ) 

YT2 = n.5^tSIN(PI*5,/6. )+STl 
DO 10CI=NFIMID,NF 
IF{S( I )-STl ) 10,10,20 
10 THITA( I ,J)=Pi/2. 

X( I , 0)=;^. 

Y( I , J)=S(I ) 

GO TO lr>o 

20 IF(S( I )-ST2)30,3C,40 
30 THITA( I , J)=5,+PI/2.-S( T )*2. 
alfa=pi/2,-thita{ r,j) 
y ( I , J)=0.5^(-( 1. -cos (ALFA) ) 

Y( T , J)=0.5*SIN(ALFA)+ST1 
GO TO Ino 

40 THITAd ,J)=-PI/3. 

X(I,J)=XT2+(S(I) -ST 2 ) ^«-COS(-THITAM) 

Y( I ,J)=YT2-(S( I ) -ST2)^^SIN!-THITAM> 
lOP CONTINUE 

1 FOPMAT(9F12.5) 

RETURN 

END 

SIDFTC FRISUP 

subroutine FRISUP( Y,THITA,V,HT,DP1,CHECK,V,NF) 
dimension Y(140,211 ,THITA(140,21) ,V(140»21) 

J = l 
1=2 

lo YP=Y(I,J) 

V( I , J) = (HT-Y{ I ,J) 

Dl=SIN(THITA(I-l,J) )/V(I-1,J) 

D2=SIN(THITA( I ,J ) )/V( I ,J) 

DY=(nl+D2)^DPl/2, 

Y( I ,J)=Y( I-l,J)4-0Y 
FR=ABS( YP-Y( I ,J) ) 

ERMAX= (HT-Y( I , J) )*CHECKV 

lF(ER-ERMAX)20,20,in 

20 CONT+NUE 

DO 3n 1=3, NF 
4r yP=y(I,J) 



V( I ,J)=(HT-Y(I ,J) )**0.5 


30 
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D1=SIN(THITA( I-2,J) )/VC I-2,J) 
D2=SIN(TH1TA(I-1»J) )/V( 

D3 = SIN(THITA( I ,J ) )/V( I ,J) 
DY=(D1+A.*D2+D3)^<-DP1/3, 

Y( I »J)=Y( 1-2 jJj+DY 
ER=ABS(YP“Y( I * J) ) 
ERMAX=(HT-Y( I » J) )*CHECKV 
lF(ER-ERMAX)3o,3n,4o 
CONTINUE 

PRINT 7, (Y( I ,J) ,I=1,N'^) 
FOi^MAT ( 10E12 .5 ) 

RETURN 

end 
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5IBFTC XYCOMP 

subroutine XYCo.>^P(X,Y,THItA,v,S,NF,NS,NFi,ms1,wFIMID,Dpi,Dvi,OP 
lDr^2,IFlRST,lTEST,H,W) 

dimension X( 140,21) ,Y( 140, 21) »THItA( i^fO, 21) ,V( 140,21 ) ,S( 140) 

FNF=NF 

PNF1=NF1 

FNS=NS 

FNS1=NS1 

I=NFIMID-1 

J = NS1 

D1=C0S(THITA(I+1,j) )/V( I+1,J) 

D2=C0S(THITA{I ,J) )/V(I,J) 

D3=COS(THITA(I-l,J) )/V{ l-l ,J) 

D= ( B.-^Dl + S ,*D2-D3 )*DP1/ 12. 

X ( I ,J)=X( I+1,U)-D 
J = NS 

Dl=-STN ( THITA( I , J-2 ) ) /V ( I , J~2 ) 

D2=-SIN( THITA( I , J-1 1 ) /V ( I , j-1 ) 

D3=-SIN(THITA( I ,J ) )/V( I ,j) 

0= (-D1+8.*D2+5.*D3)*DP1/12, 

X ( T , J)=X ( I , J-1 )-D 
I =NFIMID+1 
J = NS1 

D1=SIN(THITA{I-1,J) )/V( I-1,J) 

D2=SIN(THITA( I ,J) )/V( I ,J) 

B3=SIN(THITA( I+l,j) )/V( I+l,J) 

D= (5.*Dl+8.^fD2-D3 )^^DPl/12, 

Y(I,J)=Y(I-1,J)+D 
J = NS 

Dl=COS(THITA( I ,J-2) )/V(I,J-2) 

D2 = C0S(THITA( I ,J-1) )/V( I,J-1) 

P3=C0S(THITA( I ,J) )/V(I , J) 

D= {-Dl+8.^^D2 + 5,*D3)-i‘DPl/l2, 

S(I)=Y{I,J-1)-D 

I=NFIMID 

J=NS“2 

D1=-SIN(THITA( I,J+1) )/V{I ,J+1) 

D3=:-SIM(THITA( T,J--1) )/V(I,J-1) 

D2=-SIN(THITA( I,J) ) /VC I ,J) 

D=(5.*D1+8,*D2-D3 )*DP1/12. 

X ( I , J)=X( I ,J+1 )+D 
D1=C0S{THITA( I ,J+l) )/V( I ,J+1 ) 

D3 = C0S{ THITAC I , J-1 ) ) /V ( I , J-1 ) 

0 2=C0S(THITA( I ,J) )/V( I ,J) 

D=(5.*Dl+8,^tD2-D3}4DPl/12. 



Yd ,J)=Y( I,J + 1)+D 

I=NFlMlD-2 

J=NS 

D1=CoS(THITA( l + l,j))/vd + l,j) 
D3=C0S(THITA )/V(I-l,J) 
D2=C0 S(THITA(I,j))/V(I,j) 

D=( 5,*Dl4-8.^(-02_D3)-'^DPl/12, 

X( I »J)=X( I+l,j)-D 
Y( I » J)=o. 

I=NFIHID+2 

D1=SIN(THITA(I-1,J) )/V(I-l,J) 

D3 = SIN(THITA(I + 1,J))/VCI-*-1»J) 
D2=SIN(THITA(I ,J) )/V( I,J) 
D=(5.*D1+8,*D2-D3)*DP1/12. 

S( I )=S{ I-l)+D 
XC I » J)=n, 

J=NS 

T2=NFIMID-1 
no 15n K=3,i2 
T=NFIMID-K 
Y( I »J)=n. 
ni=l./V( 1+2, j) 

D2=1. /V( I+l, J) 

D3 = 1./V( I ,J) 

D= (D1+4,*D3+D3 )*DPl/3. 

150 X( I,J)=X( l + 2,Jl-D 
1 = 1 

no 161 J=1,NS 
FJ = J 

X ( I »J)=X ( 1,NS) 

161 Y( I ♦J)=(H+W)*( FNS-FJ)/FMS1 
1=2 

DO 155 K=l,Nsl 
J=NS-K 

IF(J-MSI) 152,151,152 

151 CONTIN'JF 

D1=-5IM(THITA{ I ,J+1) >/V(I,J+l) 
D3=-SIN(THITA( I»J-1) )/V{T,J-1) 
D2=-SIN(THITA< I,J) )/V(I,J) 

D= ( 5,*D1+8.*D2-D3 )*nPl/12. 

X ( I » J )=X ( I ,J+1 ) + 0 
D1=C0S(THITA( I ,J+1 ) )/V( I ,J+1) 
n3=C0S<THTTA( I , J-1 ) ) /V ( I , J-1 ) 
D2=C0S{THITA( I ,J) )/V(I ,J) 

D= (5,*D1+8,^D2-D3 ) »DPl/12, 



Y( I,J)=Y( I ,J+1 )+D 
GO TO 153 

152 D1=-SIN(THITA( I,J+21 )/V(I ,J+2) 
D2=-SIN(THITA( I,J+l ))/'/{! jJ+l ) 
D3=-SIN{THITA(I,J))/\/(I,J) 

D= (D1+4,*D2+D3 )-^Dpi/3o 
X(I»J)=X(I»J+2)+D 
D1 =CoS(THITA( I ,J+2) )/V( I ,J-r2) 
D2 =CoS(THITA( I ,J + 1) )/V( I ,J+1) 
D3=C0S(THITA( I jJ) )/V( I , J) 

D= (D1+4,*D2+D3 )*DP1/3. 

Y ( I > J)=Y( I , J+2 )+D 
l'^3 D1=C0S(THITA( I -1 » J ) ) / V ( I -1 , J ) 
D3=C0S(THITA( I+l , J ) ) /V ( I+l , J ) 
D2=C0S(THITA( I ,J) )/V( I ,J) 

D= (5.*D1-t- 8,*D2~D3 )*DP1/12. 
X(I»J)=(X(I-1»J)+D+X{I,J) )/2, 
D1=SIM( THITA( I-l , J) )/V( I-l , J) 
D3=SIN( THITAC I +1 » J ) ) /V ( I+l , J ) 
D2=SIN( THITA( I ,J ) )/V( I ,J) 

D= (5,*D1 + 8,-m-D2-D3 ) ;4-DP1/12. 

Y ( I > J) = ( Y( I-l» J) +D+Y ( I , J) ) /2. 
155 CONTINUE 
J = NS 

S(NFIMID)=0. 

I 3=NFIMID+3 
DO 160 1=13, NF 
nl=l./V(I-2,J) 

D2=1./V( I-1»J) 

D3=1./V( I ,J) 

D=(D1+4.*D2+D3 )*DPl/3, 

S( I )=S( I-2)+D 
160 CONTINUE 
NS2=NS-2 
DO 180 K=l,NSl 
DO 180 1=3, NF 
J=MS-K 

TF( J-NSl) 176,170>176 

176 IF ( J-NS2 ) 171 ,177, 171 

177 IF( I-NFIMID ) 171,172,171 
170 IF(I-NFIMID)173,172,173 
173 CONTINUE 

D1=-SIN (THITA( I ♦ J+1 ) J /V( I , J+i ) 
D3=-SIN(THITA( I »J-1 ) ) /V( I ,J-1 ) 
D2=-SIM (THITA( I , J) ) /V{ I ,J) 



D={ 5.*D1+8,«-D2-D3)*DP1/12. 

X( I ♦J)=X( I,J+1)+D 
D1=C0S(THITA( 1 ,J + 1) )/V( I ,J4-1) 
D3=C0S(THITA( I ,J-1) )/V( I »J-1) 
D2 = C0S(THITA( I ,jn/VC I ,J) 

D=( 5.*DH-8,«-D2-D3 )*DP1/12« 

Y( 1 »J)=Y( I ,J+1 )+D 
GO TO 172 

171 D1=-SIN(THITA( I * J+2 ) )/y( I ,J + 2) 
D2=-SIN(THITA< T »J+1 ) ) /V ( I 9 1+1) 
D3=-SIN(THITA( I »J) ) /V( I ,J) 

D=C D1+4.*D2+D3 )*DPl/3 . 

X( I 9J)=X( I ,J+2)+D 
Dl=COS(THITA( I ,J+2) )/V(I»J+2) 
D2=COS(THITA( I ,J+1 ) )/V( I ,J+l ) 
D3=C0S(THITA( I , J ) ) /V{ I ,J) 

D={ D1+4.*D2+D3 )*DPl/3, 

Y( I > J)=Y( I ,J+2)+D 

IF( ITEST.GT.O ) GO TO 180 

172 D1 = C0S(THITA( 1-2 >J) )/V(I~2»J) 
D2=C0S(THITA( I-1»J) )/V(I-1»J) 
D3=C05(THITA( I ,J) )/V( I ,J) 

D=( D1+4.*D2+D3 )*DPl/3. 

X( I 9J)=(X(I»J)+X(I-2»J)+D)/2. 
D1=SIN(THITA( I-29J) )/V(I-2,J) 
D2=SIN(THITA( I-l 9j) )/V( I-1»J) 
D3=SIN(THITA( I ,J) )/V( I ,J) 

Q= ( D14.4.*D2+D3 ) x-OPl /3, 
Y{I»J)=(Y(I-2,J)+D+Y(I ,J) )/2. 
180 CONTINUE 
RETURN 
END 
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SIBFtC input 

SUBpoUTiNE iNPUTfv twtta <- cr 

dipEnsion DYDxdonKvf 140 

f!Ml=NjFlMrD-l - ’^'^'+^’21) ,THITA{140,21) ,S(140) 

NP1=NF7MID+1 

READ <^,Nl,FMDEx 

9 EORMATI r3,5F5,2) 

1 F0RMAT(10I3) 

J=1 

w2=M1+1 

9EAD2,(DYDX(K) ,K=M2,NF) 

2 format ( 8 Flo. 5) 

read 2,(Y(I,j),i=n2,NF) 

DO 2o I=N2,NF 

20 HA ( I , j ) ( qYi^^ . j. j . 

DO 10 1=1, N1 
FI = I 
FN1=N1 

10 Y (T M ^ f h!uI^ ! TA ( M2 , J ) 4 ( ( ( F I -1 . ) /FM ) *^.FmDEX ) 

DO 30 1=1, MMl 
30 THITAII ,J)=n. 

READ 2 , ( S ( r ) , I sifiF IM rD,|'!F ) 

1 = 1 

DO 40 J=1,NS 
An THITA( I ,J)=n. 

I=NF 


DO 5n J=1,MS 
50 THITA{ I ,J)=thITAM 

print T,(THITA(I ,1) ,i = i,mF) 
PRINT 7 , { Y ( I , 1 ) , 1 = 1 ,MF) 

7 FORMAT(1oE12.5) 

return 

END 



